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New Formulas for Facilitating Osculatory 

Interpolation 

Herbert E. Salzer 

Hermite's n-point osculatory interpolation formula for equally apaced arguments at 
intervals of h, employing the function and its derivative is very much more accurate than the 
corresponding n-point Lagrangian formula and considerably more accurate than even the 
2n-point Lagrangian formula at intervals of h. Also it is specially suited for interpolation 
in many functions (e. g., Bessel, probability) that are tabulated with their derivative. To 
avoid the tremendous amount of labor in calculating the coefficients of /i and // in the forms 
that they are usually given, Hermite's formula is expressed as 



where 
and where 

being 



i i 

cti^ail(p-iy + bil{p-i),^i = ail(p — i), 

r [nl2] ]2 

[j=-[(n-l)l2] J 

[nl2] [71/2] 

n' {p-j)l n' {i-j). 

j=-[(n-l)/2] ;=-[(n-l)/2] 



The constant k(n), which may be picked arbitrarily, is here chosen to make a^ and &,• integers. 
The exact values of a^ and bi are given for n = 2{l)ll, i=—[(n—l)l2] to [n/2] so that this 
formula can be applied exactly for any polynomial up to the 21 St degree. A schedule gives 
approximate upper bounds for the coefficients of /(2n)(^j/j,2n^^2nj(^^) in j^ (p)^ 

When a function /(a;) and its first derivative are known at n points Xf, i=l, 2, . . . ^ n, a 
highly accurate interpolation formula due to Hermite is given by 

f{x) = ±{Li^\x)}^l-2L^-^\x,)(x-xd]f(Xi)+±{Lt\xmx-x,)n^^^ (1) 

where 

n n 

Ll'^\x)= ll'{x — Xj)l\l\Xi—Xi), j=i is absent from n', (2) 

j=i j=i 

and 

^2n(^) = /^'"^(?)|n(x-x^)| /{2n)\, least Xi<^< greatest Xi. (3) 

Thus (1) is exact whenever /(ic) is a polynomial of degree less than or equal to 2n—l. 

When the points Xt are equally spaced at intervals of A, it is customary to alter the notation 
in Xi, letting i run from — [(n— 1)/2] to [n/2] instead of 1 to n, where [m] denotes the largest 
integer not exceeding m. Then it is convenient to choose a variable p given by x=XQ+ph and 
to let Xi=Xo+ih. Also, f(x)=f(xo+ph)=fp=f, J{:Xi)=Ji, SindfiXi)=fi. Then (1), for /(a;) 
considered as a function of p, is expressible as 

[n/2] [nf2] 

f{xo+ph)= s {z^>(p)}^{i-2z^>'(i) {p-i)}fi+ s [Lnvmp-i)hf:+R2n(p\ 

i=-[(n-l)/2] i=-[(n-l)l2] 

(4) 
where now 

N21 [nl2] 

Lt\v)= n' (p-iV n' {i-j), (5) 

;=-[(n-l)/2] i=-[(n-l)/2] 

and 

( {nl2] \ 2 

R2n{v)-S''''\iW\ n {p-j)\li2n)\, a;_u„_,),,,<j<xi„/2i. (6) 

(j=-r(7i-l)/2] ) 
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There are many advantages in the use of (1) or (4) over the ordinary Lagrangian interpola- 
tion formula given (for equal spacing) by 

}{x,+vh)= S Ll-\v)fi + Rniph (7) 

i=-[{n-l)l2] 

where 

[n/2] 

Rn(p)=f''''m^ n {p-j)/nl, X_[(n-l)/21<?<aJ[n/2]. (8) 

j=-[(n-l)/2] 

N2] 

Thus letting 11 (p—j) be denoted by L^'^^(p), and recalling that for reasonably small h 

; = -[(n-l)/2] 

we have approximately 

f'^\^)h'^^A% (9) 

where A'^J is the approximate mth difference of the tabulated f{;x) , the remainder term for (7) 
is of the order of A'^JL^''\p)ln\, whereas that for (4) is of the order of A^'^j{L^'^\'p)Yl{2n)\. 
Apart from the fact that A^^^j is usually very much smaller than A^, the factor [L^""^ {p)Y j {2n)\ 
is equal to the square of L^''\p)ln\ (where the latter is usually very small and less than unity 
so that its square is ever so much smaller) multiplied by the very small quantity 7i!/(n+l) 
. . . {2n). The user can appreciate the improvement by comparing the approximate upper 
bounds for the multiplier oi J^^''\^)h^'' in R2n{p) of (4) which are tabulated in the schedule at 
the end, with the approximate upper bounds for the multiplier of /^''^(J)/i'* in Rn{p) of (7), 
which are tabulated in [3, p. xvi]^ and from which this present schedule was calculated. Thus 
it will be apparent that (4) is a very much more accurate formula than (7). Of course, we are 
comparing (4), a confluent form of a 27i-point formula, with (7), which is only an n-point formula. 

But it is important to note than even if (4) is compared with formula (7) taken for 2n 
points at intervals of A, instead of n points, the remainder term would differ from that in (4) 
(apart from a different J in/^^"^($)) by the presence of the factor L^'^^^p) instead of the factor 
|2;<«)(^)}2^ which has a very much smaller upper bound than the Z^^"^(7)), showing that wher- 
ever it is possible to be used, the 7i-point Hermite osculating interpolation formula is much 
to be preferred, regarding accuracy, to a 27?-point Lagrangian interpolation formula at 
the same interval h. This last statement becomes intuitively plausible when the osculating 
interpolation formula for n points at intervals of h is regarded as a confluent form of a 
2n-point Lagrangian formula whose 2n points lie within a range of nh so that the ^'average inter- 
vaP' between those 2n points is only half the interval of h for the usual 2?i-point Lagrangian 
formula. Thus the upper bound for the remainder term of the n-point osculating formula 
would be expected to be of the order of (l/2^^)th of the upper bound for the remainder term of 
the 271-point Lagrangian formula; actual estimates show it to be even considerably smaller. 
For example, the 2-, 3-, 4-, and 5-point osculating formulas have error terms whose upper 
bounds are around l/16th, 1/110, 1/640, and 1/3000 of the respective upper bounds of the error 
terms in the 4-, 6-, 8-, and 10-point Lagrangian formulas. 

A second important advantage in (4) is that it is suited for use with very many tables 
where the derivative of the function is tabulated alongside of the function itself. For example, 
it is useful in tables of Bessel functions of the first and second kind [4, 5, 6] which give Ji{x) = 
—Jq(x)j Yi{x) = — Yq{x), and probabilit}^ functions [7], and in numerous tables of more ele- 
mentary functions and their integrals, such as tables of sine, cosine, or exponential integrals 
[8, 9], where the derivative is very easy to obtain. 

However, the use of (1) or (4) in the form usually presented [1, 2] requires a considerable 
amount of computational labor which mounts considerably as the number of points Xi 
increases. It is the purpose of this present article to provide a means of using (4) with a small 
fraction of the labor involved in the direct calculation of the coefficients of /^ and/-. The idea 

* Figures in brackets indicate the literature references at the end of this paper. 
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behind this method goes back to a scheme first used b}^ W. J. Taylor for (•alciilatiiig Tjagraiigiari 
interpolation coefficients for functions tabulated at real equidistant arguments [10], and which 
was generalized by the present writer for functions tabulated at nonequidistant arguments 
[11], and also for complex arguments whether in Cartesian [12] or polar form [13, 14], and 
finally even for functions that are interpolable by expressions that are not transformable into 
polynomials [15]. Recently, the writer in lool^ing for some way to reduce the amount of work 
in using (4), observed that Taylor's idea could be extended also to the calculation of osculating 
interpolation coefficients. In place of extensive tables of the (2/^— l)th degree polynomial 
coefficients of /^ and /• in (4) , one requires for each separate n only some fixed quantities ai 
and bi, which are exact integers and are tabulated below. To see this, one merely expresses 
(4) as 

[nl2\ (/ A2 9 r ("> V'}*^ /4 2\ 42 \ 

/(x.+ M)=(X"-Wl-,._2„J(i^,-4^p)/,+g^A/;| + «.W, (10) 
where ^ 

/ [n/2\ 

Ai^i n' {i-j), (11) 

/i=-[(n-i)/2] 

Now the right member of (4) or (10) without the R2n{p) gives the expression for a (2n—l)th. 
degree polynomial, which, with its derivative, assumes preassigned values of /^ and/- at x=Xu 
and moreover that polynomial is uniqyiAy determined by the ft aridj[. For proof of uniqueness 
see [1, p. 85-86], where T. Fort gives a demonstration of the unique existence of a more general 
osculating formula. His proof is practically complete save for the explicit indication that 
the mode of representation of any (mn— l)th degree polynomial which is given at the bottom 
of page 85 is always possible (which is fairly obvious). Now we make use of this uniqueness 
of representation by putting /(.t) = 1 into (10), so that both/- and Ronip) are zero,/j=l, and 
we get 

{L''ip)V-—u^, X ,/ oru,>/.^^.2^ • (12) 



1^1 / A^ 2Lt''{i)A^ 

•[(ri-l)/2lV' 

Thus from (10) and (12), 



j=_[(„-l)/21 



Kip-if p-i ) 



i=_[tr-l)/2]\(2?-^)' (2?-W 



where at and bi are given by 



a, = k{n)A\, (14) 

b, = k{n){-2Lr(i)A\], (15) 



and where k(jib) is any suitably chosen constant of proportionality that depends only upon n\ 
In the present case the k(ji) was chosen as to give exact integral values for a^ and bi instead 
of rational fractional values. 

It is simplest to think of the approximation to j{x) in the concise form 

f^^S-ifi+lMl, (16) 

where 

at^ai/(p-ir + bi/(p-i), (17) 

and 

p.^a,l(p-i). (18) 



a The dependence upon n of Ai, as well as of a» and bi given below, is not indicated, so as to avoid cumbersome notation. 
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In using (16), (17), and (18) with a desk calculator, it is easiest to first divide ai by jp—i to 

get Pi, which is next both multiphed by hf/ and increased by hi. The latter, or Pi+hi, is again 

divided by {p—i) to give a^, from which one obtains both a^/j and Sa^ and fin ally 2(a:i/i+/3iA//)/Sai. 

The computation of the quantities at and hi was quite straightforward. Since Ai= 

(_l)[w/2]-u ^ j /(^_1)!^ instead of A^, the proportional quantities ( . , ./ -.n /m ) 

\i+[{n— I) 12]/ 1 \^+ [(n— 1)/2J/ 

were calculated. Then they were multiplied by the — 2Z/'^^^(i), which were calculated by 
differentiating the explicit polynomial expressions Li^'^^p) and then setting p=i. All frac- 

tionsin— 21 . , ,^ x, t) ii^'^^'('i) were cleared by multiplication of these quantities, as well as 

\^+[(n~l)/2]/ 

the ( . ,, X , -, I > hy some suitable integer, for each n, to yield the exact integral values 

\^+[(n~■l)/2]/ * 

for Gi and 6^, which are tabulated below. The at and 6^ were checked by both recomputation 
and by use in an example for every n where the answer was known exactly and where the com- 
putation by (13) (or the equivalent (16), (17), and (18)) doing the work in decimal form to 
avoid too much labor, gave agreement to 10 significant figures. 

The schedule giving the approximate upper bounds for the coefficients of f^^"^^ (^)h^^ ^ A^^f 
in the error term B>2n{p), (see (4) with (6), or (13)), namely, the quantities [L^'^\p)YI{2n)\, 
was calculated from the approximate upper bounds for L^^\p)/nl given in a schedule in [3, p. 
xvi], by squaring the entries in the latter and multiplying by nl/(n-{-l) . . . {2n). Thus, 
since the L^'^\p)/nl was tabulated only approximately, in some cases to only one significant 
figure, some of the upper bounds given here for {L^''^(p)}^/(2n)l are not at all precise. For 
example, if we take a rounded 0.01 and square it to obtain 0.0001, the true value of that square 
may be only one-fourth as large or over twice as large, depending upon whether the 0.01 was 
rounded from 0.0051 or 0.0149. Hence in the schedule below, in some cases for the larger 
number of points, only the order of magnitude of an upper bound for the {L^'^^(p)}^/(2n)l is 
indicated. But more precise determination will hardly be needed there due to the extreme 
accuracy that is surely indicated even within the range of the uncertainty. 
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Table of a,- and hi 



^^ 





Two-point 






Seven-point 






Ten-point 




ao 


1 


6o 


2 


a-3 


10 


6-3 


49 


a-i 


1260 


6-4 


7129 


«i 


1 


h 


-2 


a-2 


360 


6-2 


924 


a-3 

1 


1 02060 


6-3 


3 50649 








a-i 


2250 


6-1 


2625 ' 


a-2 


16 32960 


6-2 


35 69184 




Three-point 




ao 


4000 


60 





a-i 


88 90560 


6-1 


109 65024 


a-i 


1 


h-x 


3 


Ol 


2250 


61 


-2625 


ao 


200 03760 


60 


80 01504 


flo 


4 


ho 





ai 


360 


62 


-924 


ai 


200 03760 


61 


-80 01504 


ai 


1 


hi 


-3 


«3 


10 


63 


-49 


a-i 


88 90560 


62 


-109 65024 














ai 


16 32960 


63 


— 35 69184 




Four-point 






Eight-point 




a^ 


1 02060 


64 


-3 50649 


a-i 


3 
27 


h-i 

ho 


11 
27 


a-3 


70 


6-3 


363 


as 


1260 


65 


-7129 


ao 


a-2 


3430 


6-2 


9947 








ai 


27 


h 


-27 


a-i 


30870 


6-1 


48363 




Eleven-point 




a2 


3 


h2 


-11 


ao 


85750 


^o 


42875 


a-5 


1260 


6-5 


7381 








ai 


85750 


61 


— 42875 


a-4 


1 26000 


6-4 


4 60900 




Five-point 




a2 


30870 


62 


-48363 


a-3 


25 51500 


6-3 


62 14725 


a-2 


6 


6-2 


25 


az 


3430 


63 


-9947 


a-2 


181 44000 


6-2 


275 61600 


a-i 


96 
216 


6-1 
ho 


160 



04 


70 


64 


-363 


a-i 


555 66000 


6-1 


407 48400 


ao 








ao 


800 15040 


60 





ai 


96 


h 


-160 




Nine-point 




ai 


555 66000 


61 


-407 48400 


a2 


6 


h2 


-25 


a-4 


140 


6-4 


761 


a2 


181 44000 


62 


-275 61600 








a-3 


8960 


6-3 


28544 


a3 


25 51500 


63 


-62 14725 




Six-point 




a-2 


1 09760 


6-2 


2 08544 


ax 


1 26000 


64 


-4 60900 


a_2 


30 


h..2 


137 


a-i 


4 39040 


6-1 


3 95136 


a5 


1260 


6. 


-7381 


a-i 


750 


6-1 


1625 


ao 


6 86000 


60 











Go 


3000 


6o 


2000 


ai 


4 39040 


61 


-3 95136 








ai 


3000 


6, 


-2000 


a2 


1 09760 


62 


-2 08544 








a2 


750 


62 


-1625 


as 


8960 


63 


- 28544 








dz 


30 


63 


-137 


a4 


140 


64 


-761 









Schedule of approximate upper hounds for {L^">(p) \^/{2n)\ 
(Figures in parentheses indicate the number of zeros between the decimal point and the first significant digit.) 



Range of p__ 
Two-point--- 
Four-point__ 

Six-point 

Eight-point-. 
Ten-point 

Range of /?__ 

Three-point _ 
Five-point- -- 
Seven-point _ 
Nine-point- _ 
Eleven-point 



0<P<1 



-1<P<0 
1<P<2 



-2<p<-l 
2<p< 3 



-3<p<-2 
3<p< 4 



-4<p<-3 
4<p< 5 



.0026 

. (5)82 
. (7)26 
. (10)94 
. (12)34 



(4)25 
(7)55 
(9)15 
(12)49 



(6)62 
(9)85 
(11)18 



(7)20 
(10)18 



(9)7^ 



o<bl<i 



i<bl<2 



2<\p\<S 



3<ip|<4 



4<bl<5 



. 00021 

.(6)57 
.(8)18 
. (11)60 
. (13)6 



(5)38 
(8)62 
(10)15 

(13)6 



. (6)11 
. (9)12 
. (12)2 



(8)40 
(11)6 



(9)1 



Washington, May 12, 1953. 
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